c  Tzong-Yow Hwu
c Fortran routines for use in p3d reader for fortran unformatted data files
      subroutine foropn(lun, iname, length, opstat)
	integer lun, iname(1800), length, opstat 
c       lun: the logical unit number for the opened file
c       iname: the ascii code of the file name
c       length: the length of the iname integer array
c       opstat: an error indicator, a value of 0 indicates an error on open
c             and other values O.K.
c
	character*180 sname
	integer done, notdone
c
	parameter(done = 1, notdone = 0)
	opstat = notdone 
c
c Mapping the file name to sname
	do istep=1,length
          sname(istep:istep)=char(iname(istep))
        end do
        open(unit=lun, err=100, file=sname(1:length), status='old',
     &       access='sequential', form='unformatted')
	opstat = done
100     return
c
      end
c
c
      subroutine forcls(lun)
	integer lun
c
	close(unit=lun)
	return
      end
c
c
c  For reading ngrid in fortran unformatted file
      subroutine rngrid(lun, ngrid, opstat) 
	integer lun, ngrid, opstat 
c       lun:     logical unit number connected to the opened xyz file
c       ngrid:   ngrid to be read
c       opstat:  error indicator
c
	integer*4 status
	integer done, notdone, eofile
	parameter(done = 1, notdone = 0, eofile=-1)
c
	opstat = notdone 
	read(unit=lun, err=200, iostat=status) ngrid
c
	opstat = done 
200     if (status .LT. 0) then
	  opstat = eofile
        end if
        return
      end
c
c
c  to use this routine, indims must be allocated as an int array of size
c  ngrid*ndim where ndim is the number of dimensions for xyz and solution
c  file, and number of dimensions plus one(for nvar) for function file
c  For reading indims in fortran unformatted file
      subroutine rddims(lun, ndim, ngrid, indims, opstat) 
	integer lun, ndim, ngrid, opstat 
	integer indims(ndim*ngrid)
c       lun:     logical unit number connected to the opened xyz file
c       ndim:    number of dimensions
c       ngrid:   ngrid to be read
c       indims:  dimension array
c       opstat:  error indicator
c
	integer*4 status
	integer done, notdone, eofile
	parameter(done = 1, notdone = 0, eofile=-1)
	integer i, j, k, n
c
	opstat = notdone 
	read(unit=lun, err=300, iostat=status) 
     &       ((indims(n*ndim+i), i = 1, ndim), n = 0, ngrid - 1)
c
	opstat = done 
300     if (status .LT. 0) then
	  opstat = eofile
        endif
        return
      end
c
c
c  For reading grid values in fortran unformatted file
c  Reading a grid values of a single grid
      subroutine rdgrid
     &(lun, ndim, indims, isiblk, iblank, iperm, gdvals, size, opstat) 
	integer lun, ndim, indims(ndim), isiblk, iperm, opstat
	integer*4 size
	integer iblank(size)
	real gdvals(size*ndim) 
c       lun:     logical unit number connected to the opened xyz file
c       ndim:    number of dimensions
c       indims:  indirect dimension array
c       isiblk:  is there a iblank in the grid file
c       iblank:  iblank array
c       iperm:   the arrangement of the grid value in whole or plane
c       gdvals:  grid values
c       size:    the size of one component of the grid
c       opstat:  error indicator
c
	integer done, notdone, whole, yes
	parameter(done = 1, notdone = 0, whole=0, yes=1)
	integer i, j, k, n
c
	opstat = notdone 
c    
        if (iperm .EQ. whole) then
c       grid data is arranged in whole
	  if (isiblk .EQ. yes) then
c       grid file contains iblank values
            if (ndim .EQ. 1) then
              read(unit=lun, err=400) 
     &          (gdvals(i), i=1, indims(1)), 
     &          (iblank(i), i=1, indims(1))
	    else if (ndim .EQ. 2) then
	      read(unit=lun, err=400)
     &        (((gdvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                             j=1,indims(2)),n=0,ndim-1),
     &        ((iblank(i*indims(2)+j), i=0,indims(1)-1), j=1,indims(2))
	    else if (ndim .EQ. 3) then
	      read(unit=lun, err=400)
     &        ((((gdvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &        i=0,indims(1)-1),j=0,indims(2)-1),k=1,indims(3)),
     &        n=0,ndim-1),
     &        (((iblank(i*indims(2)*indims(3)+j*indims(3)+k),
     &        i=0,indims(1)-1),j=0,indims(2)-1),k=1,indims(3))
	    else
	      go to 400
            end if
c
	  else
c
c       grid file contains no iblank values
c
            if (ndim .EQ. 1) then
              read(unit=lun, err=400) 
     &          (gdvals(i), i=1, indims(1)) 
	    else if (ndim .EQ. 2) then
	      read(unit=lun, err=400)
     &        (((gdvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                             j=1,indims(2)),n=0,ndim-1)
	    else if (ndim .EQ. 3) then
	      read(unit=lun, err=400)
     &        ((((gdvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &        i=0,indims(1)-1),j=0,indims(2)-1),k=1,indims(3)),
     &        n=0,ndim-1)
	    else
	      go to 400
            end if
	  endif
c
	else
c       grid data is arranged in plane
c
	  if (isiblk .EQ. yes) then
c       grid file contains iblank values
            if (ndim .EQ. 1) then
              read(unit=lun, err=400) 
     &          (gdvals(i), i=1, indims(1)), 
     &          (iblank(i), i=1, indims(1))
	    else if (ndim .EQ. 2) then
	      read(unit=lun, err=400)
     &        (((gdvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                             j=1,indims(2)),n=0,ndim-1),
     &        ((iblank(i*indims(2)+j), i=0,indims(1)-1), j=1,indims(2))
	    else if (ndim .EQ. 3) then
	      do k = 1, indims(3)
	        read(unit=lun, err=400)
     &          (((gdvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &          i=0,indims(1)-1),j=0,indims(2)-1), n=0,ndim-1),
     &          ((iblank(i*indims(2)*indims(3)+j*indims(3)+k),
     &          i=0,indims(1)-1),j=0,indims(2)-1)
	      end do
	    else
	      go to 400
            end if
c
	  else
c
c       grid file contains no iblank values
c
            if (ndim .EQ. 1) then
              read(unit=lun, err=400) 
     &          (gdvals(i), i=1, indims(1)) 
	    else if (ndim .EQ. 2) then
	      read(unit=lun, err=400)
     &        (((gdvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                             j=1,indims(2)),n=0,ndim-1)
	    else if (ndim .EQ. 3) then
	      do k = 1, indims(3)
	        read(unit=lun, err=400)
     &          (((gdvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &          i=0,indims(1)-1),j=0,indims(2)-1), n=0,ndim-1)
	      end do
	    else
	      go to 400
            end if
	  endif
c
	endif
c
	opstat = done 
400     return
      end
c
c
c  For reading solution values in fortran unformatted file
c  Reading solution values of a single grid
      subroutine rdsolu
     &(lun, ndim, indims, iperm, slvals, size, opstat) 
	integer lun, ndim, indims(ndim), iperm, opstat
	integer*4 size
	real slvals(size*(ndim+2)) 
c             ndime+2 since includes density and pressure
c       lun:     logical unit number connected to the opened Q file
c       ndim:    number of dimensions
c       indims:  indirect dimension array
c       iperm:   the arrangement of the solution values in whole or plane
c       slvals:  density, pressure, and solution values
c       size:    the size of one component of the solution values 
c       opstat:  error indicator
c
	integer done, notdone, whole
	parameter(done = 1, notdone = 0, whole=0)
	integer i, j, k, n
c
	opstat = notdone 
c    
        if (ndim .EQ. 1) then
          read(unit=lun, err=500) 
     &      ((slvals(n*size+i), i=1, indims(1)), n = 0, 2) 
	else if (ndim .EQ. 2) then
	  read(unit=lun, err=500)
     &    (((slvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                         j=1,indims(2)),n=0,3)
	else if (ndim .EQ. 3) then
          if (iperm .EQ. whole) then
c           solution data is arranged in whole
	    read(unit=lun, err=500)
     &      ((((slvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &      i=0,indims(1)-1),j=0,indims(2)-1),k=1,indims(3)),
     &      n=0,4)
	  else
c           solution data is arranged in plane
	    do k = 1, indims(3)
	      read(unit=lun, err=500)
     &        (((slvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &        i=0,indims(1)-1),j=0,indims(2)-1), n=0,4)
	    end do
	  end if
	else
	  go to 500
        end if
c
	opstat = done 
500     return
      end
c
c
c     for use to read time data from solution q files
      subroutine rdtime(lun, time, opstat)
	integer lun, opstat
	real time
c       lun:     logical unit number connected to the opened Q file
c       time:    the time of the dataset
c       opstat:  error indicator
c
c       useless data to be discarded
        real fsmach, alpha, re 
	integer done, notdone, whole
	parameter(done = 1, notdone = 0, whole=0)
c
	opstat = notdone 
	read(unit=lun, err=600) fsmach, alpha, re, time
        opstat = done	
600     return
      end
c
c
c  For reading function values in fortran unformatted file
c  Reading function values of a single grid
      subroutine rdfunc
     &(lun, ndim, indims, nvar, iperm, fnvals, size, opstat) 
	integer lun, ndim, indims(ndim), nvar, iperm, opstat
	integer*4 size
	real fnvals(size*nvar) 
c       lun:     logical unit number connected to the opened function file
c       ndim:    number of dimensions
c       indims:  indirect dimension array
c       nvar:    rank: value of 1 means scalar, more means vector 
c       iperm:   the arrangement of the solution values in whole or plane
c       funcvals:  function values of n variables 
c       size:    the size of the function values 
c       opstat:  error indicator
c
	integer done, notdone, whole
	parameter(done = 1, notdone = 0, whole=0)
	integer i, j, k, n
c
	opstat = notdone 
c    
        if (ndim .EQ. 1) then
          read(unit=lun, err=700) 
     &      ((fnvals(n*size+i), i=1, indims(1)), n = 0, nvar-1) 
	else if (ndim .EQ. 2) then
	  read(unit=lun, err=700)
     &    (((fnvals(n*size+i*indims(2)+j), i=0,indims(1) - 1),
     &                         j=1,indims(2)),n=0,nvar-1)
	else if (ndim .EQ. 3) then
          if (iperm .EQ. whole) then
c           solution data is arranged in whole
	    read(unit=lun, err=700)
     &      ((((fnvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &      i=0,indims(1)-1),j=0,indims(2)-1),k=1,indims(3)),
     &      n=0,nvar-1)
	  else
c           solution data is arranged in plane
	    do k = 1, indims(3)
	      read(unit=lun, err=700)
     &        (((fnvals(n*size+i*indims(2)*indims(3)+j*indims(3)+k),
     &        i=0,indims(1)-1),j=0,indims(2)-1), n=0,nvar-1)
	    end do
	  end if
	else
	  go to 700
        end if
c
	opstat = done 
700     return
      end
